function  [OB,OR,ORU,Orders,Est] = UCombMAX(D,PL,PR,alpha,B,Orders)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% File:               UCombMAX.m
%
% Author:             Tasos Kalandrakis
%
% Description:        Computes count estimator of order \pi and builds
% confidence sets using test procedure of Romano et al (2014) with the MAX,
% statistic and moment inequalities (M_B), (M_R), or (M_R)+(M_U).
% Uses square sample.
%
% Dependencies:       Calls 'ModMAXTestTwoStep.m' and 'RUMoment.m'
%
% Created:            Dec - 08 - 2021
%
% Last Modified:      March - 1 - 2022
%
% Language:           MATLAB
%
% Related Reference: Tasos Kalandrakis. "One-dimensional scaling without
%                    without apologies," forthcoming, Journal of Politics
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sel = (sum(isnan(D),2)==0);
D = D(sel,:);
[N,P] = size(D);

if nargin<5
    B=499;
end
if nargin<6
    Orders=perms(1:P);
    Norders=length(Orders(:,1));
    MarkO=zeros(Norders,1);
    for o=1:Norders
        for p=1:P
            if Orders(o,p)==PR
                PRloc=p;
            elseif Orders(o,p)==PL
                PLloc=p;
            end
        end
        if PRloc<PLloc
            MarkO(o,1)=1;
        end
    end
    Orders=Orders(MarkO==0,:);
end
Norders=length(Orders(:,1));

Triplets=P*(P-1)*(P-2);
Appears=zeros(Norders,Triplets);
S1=zeros(N,Triplets);
S2=zeros(N,Triplets);
BM1=zeros(N,Triplets);
BM2=zeros(N,Triplets);
RM=zeros(N,Triplets);
col=1;
Ranks=repmat([1:P],Norders,1);
for r1=1:P
    for r2=setdiff(1:P,r1)
        for r3=setdiff(1:P,[r1 r2])
            Appears(:,col)=sign(sum(Ranks.*(Orders==r2),2)-sum(Ranks.*(Orders==r1),2))+sign(sum(Ranks.*(Orders==r3),2)-sum(Ranks.*(Orders==r2),2))==2;
            S1(:,col)=sign(D(:,r1)-D(:,r2));
            S2(:,col)=sign(D(:,r2)-D(:,r3));
            BM1(:,col)=(D(:,r1)<min([D(:,r2) D(:,r3)],[],2))-...
                (D(:,r2)<min([D(:,r1) D(:,r3)],[],2));
            BM2(:,col)=(D(:,r3)<min([D(:,r1) D(:,r2)],[],2))-...
                (D(:,r2)<min([D(:,r1) D(:,r3)],[],2));
            RM(:,col)=(D(:,r1)<min([D(:,r2) D(:,r3)],[],2))+...
                (D(:,r3)<min([D(:,r1) D(:,r2)],[],2))-...
                2*(D(:,r2)<min([D(:,r1) D(:,r3)],[],2));
            col=col+1;
        end
    end
end

SD=(S1-S2==2);
Triscore=sum(SD);
criterionQ=Appears*Triscore';

RM=-RM;
BM=-[BM1 BM2];

Appears2=[Appears Appears];
trips=Triplets/6;

%storage of statistics
sB = zeros(Norders,1);
sR = zeros(Norders,1);
sRU = zeros(Norders,P);

rejectB=zeros(Norders,1);
rejectR=zeros(Norders,1);
rejectRU=zeros(Norders,P);

critVB=zeros(Norders,1);
critVR=zeros(Norders,1);
critVRU=zeros(Norders,P);

muNullB=zeros(Norders,trips*2);
muNullR=zeros(Norders,trips);
muNullRU=cell(Norders,P);

parfor i=1:Norders
    %B
    BMi=BM(:,Appears2(i,:)==1);
    [rejectB(i,:),sB(i,:),critVB(i,:),muNullB(i,:)] = ModMAXTestTwoStep(BMi,alpha,B);
    
    %R
    RMi=RM(:,Appears(i,:)==1);
    [rejectR(i,:),sR(i,:),critVR(i,:),muNullR(i,:)] = ModMAXTestTwoStep(RMi,alpha,B);
    
    %RU
    for j=1:P
        RUMij = RUMoment(D,Orders(i,:),j);
        [rejectRU(i,j),sRU(i,j),critVRU(i,j),muNullRU{i,j}] = ModMAXTestTwoStep(-RUMij,alpha,B);  
    end

end


IndOrders=1:Norders;

Est.N=N;
Est.P=P;
Est.Orders=Orders;

Est.criterionQ=criterionQ;
Est.Qsc=min(criterionQ);
Est.Qorder=Orders(criterionQ==Est.Qsc,:);
Est.Qind=IndOrders(criterionQ==Est.Qsc);

OB.reject=rejectB;
OB.critVal=critVB;
OB.muNull=muNullB;
OB.test=sB;

OR.reject=rejectR;
OR.critVal=critVR;
OR.muNull=muNullR;
OR.test=sR;

ORU.reject=rejectRU;
ORU.critVal=critVRU;
ORU.muNull=muNullRU;
ORU.test=sRU;

end